## ---- simulation-functions ----

trunctedProp <- function(n, mean, sd = .25) {
  require(truncnorm)
  return(rtruncnorm(n=n, a=0, b=1, mean=mean, sd = sd))
}

randomSocialNetwork <- function(n_communities, g_size, min_comm_size,
                                within_p, between_p) {
  
  if(n_communities>20) stop("Up to 20 communities")
  
  library(igraph)
  library(openssl)
  
  g_list <- list()
  default_comm_size <- g_size / n_communities
  comm_sizes <- rep(default_comm_size, n_communities)
  
  # Randomise community size
  isOdd <- function(x) {x %% 2 != 0}
  for (i in 1:n_communities) {
    if (isOdd(i)) {
      if (i == n_communities) next
      this_change <- 
        sample(-abs(default_comm_size-min_comm_size):abs(default_comm_size-min_comm_size), 
               1)
      comm_sizes[i] <- comm_sizes[i] + this_change
      last_change <- -this_change
    } else {
      comm_sizes[i] <- comm_sizes[i] + last_change
    }
  }
  
  for (i in 1:n_communities) {
    # print(paste0(i, " ", comm_sizes[i]))
    g_list[[i]] <- random.graph.game(comm_sizes[i], within_p)
    V(g_list[[i]])$community <- LETTERS[i]
    V(g_list[[i]])$name <- paste0(LETTERS[i], 1:vcount(g_list[[i]]))
  }
  
  g <- do.call(disjoint_union, g_list)
  v_communities <- V(g)$community
  g_mat <- as_adj(g)
  from_is <- c(1,cumsum(comm_sizes)+1)
  from_is <- from_is[1:n_communities]
  to_is <- cumsum(comm_sizes)
  is <- mapply(seq, from_is, to_is, by = 1)
  for (i in 1:n_communities) {
    g_mat[is[[i]], unlist(is[which(1:n_communities != i)])] <- 
      rbinom(n=length(as.numeric(g_mat[is[[i]], 
                                       unlist(is[which(1:n_communities != i)])])), 
             size=1, prob=between_p)
  }
  g <- graph.adjacency(g_mat, mode = 'undirected')
  V(g)$community <- v_communities
  return(g)
}

activationProb <- function(political_distress, 
                           internal_drive, relational_drive, universal_drive,
                           neighborhood, i) {
  p <- 
    1/2 * (political_distress * internal_drive) + 
    1/2 * (political_distress * (relational_drive ^ (1/neighborhood)))
  return(p) 
}

addSCNet <- function(g, prob = .3, mean = 1, sd = .25) {
  nodes <- V(g)$name[V(g)$internal_drive == 1]
  comb_nodes <- combn(nodes, 2)
  comb_nodes <- comb_nodes[,as.logical(rbinom(ncol(comb_nodes), size = 1, prob))]
  g <- g + edges(as.character(comb_nodes))
  g <- simplify(g)
  edge_attr(g, 'weight', get.edge.ids(g, as.character(comb_nodes))) <-
    trunctedProp(ncol(comb_nodes), mean=1, sd = .1)
  return(g)
}

getWeightedNeigh <- function(neigh, g) {
  len <- length(neigh)
  if (len == 1) return(0)
  ego <- neigh$name[1]
  neighbors <- neigh$name[2:len]
  active <- neigh$active[2:len]
  ids <- get.edge.ids(g, as.vector(rbind(rep(ego, len-1), neighbors)))
  weight <- edge_attr(g, name = 'weight', index = ids)
  return(weighted.mean(active, weight))
}

updateNetwork <- function(g, i, stay_active = TRUE, reset_each = NULL) {
  neigh_list <- neighborhood(g, order = 1)
  neigh <- unlist(lapply(neigh_list, getWeightedNeigh, g))
  neigh[is.na(neigh)] <- 0
  prob <- activationProb(V(g)$political_distress,
                         V(g)$internal_drive, 
                         V(g)$relational_drive, 
                         V(g)$universal_drive,
                         neigh, i)
  active <- as.logical(rbinom(vcount(g), size = 1, p = prob))
  if (stay_active == TRUE) {
    active[V(g)$active == TRUE] <- TRUE
  } else {
    active[V(g)$active == TRUE & V(g)$past < reset_each] <- TRUE
  }
  past <- ifelse(active == FALSE, 0, V(g)$past + active)
  return(list('active' = active, 'past' = past))
}

doNetSim <- function(i, g, path = NULL, print = FALSE) {
  if(i > 1) {
    res <- updateNetwork(g, i, stay_active = T)
    V(g)$active <- res[['active']]
    V(g)$past <- res[['past']]
  }
  if (print == TRUE) {
    jpeg(file = paste0(path, sprintf('%03d.jpg', i)),
         width = 2500, height = 2500)
    par(mar=c(0,0,4,0))
    plot(g, vertex.label = NA, vertex.size = 4, layout = g_lo,
         vertex.color = ifelse(V(g)$active, "black", "white"))
    title(main = paste0("Time: ", i, 
                        "     Activation: ", 
                        round(sum(V(g)$active) / vcount(g) *100, 2), 
                        "%"),
          cex.main = 6)
    dev.off()
  }
  return(list(graph = g, activation = sum(V(g)$active) / vcount(g)))
}